pacman::p_load('tidyverse',
'sf',
'here',
'janitor',
'tmap',
'kableExtra',
'patchwork',
'ggthemes',
'colorRamps',
'stars',
'terra',
'leaflet',
'leaflet.extras')Leaflet Exploratory
Load packages
Read Rast (calfire)
rast_path <- here::here('data','ds1327', 'tiff', 'ds1327.tif')
gdb_path <- here::here('data','ds1327', 'ds1327.gdb')
habitat <- rast(gdb_path)veg_type <- habitat
activeCat(veg_type, layer = 1) <- 'LIFEFORM'
veg_type <- project(veg_type, "EPSG:4326") %>%
aggregate(fact = 20)
|---------|---------|---------|---------|
=========================================
|---------|---------|---------|---------|
=========================================
Read GAP GDB
gap <- st_read(here::here('data',
'PADUS4_1_State_CA_GDB_KMZ',
'PADUS4_1_StateCA.gdb'),
quiet = TRUE) %>%
clean_names() Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
automatically selected the first layer in a data source containing more than
one.
Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : GDAL
Message 1: organizePolygons() received a polygon with more than 100 parts. The
processing may be really slow. You can skip the processing by setting
METHOD=SKIP, or only make it analyze counter-clock wise parts by setting
METHOD=ONLY_CCW if you can assume that the outline of holes is counter-clock
wise defined
sf_use_s2(FALSE)Spherical geometry (s2) switched off
gap_clean <- gap %>%
select(own_type, gap_sts, mang_type, Shape) %>%
st_cast("MULTIPOLYGON") %>% # Convert curves FIRST
st_make_valid() %>%
st_transform(4326) %>%
st_simplify(dTolerance = 0.01) %>%
group_by(gap_sts) %>%
summarise(geometry = st_union(Shape), .groups = 'drop')Warning in st_simplify.sfc(st_geometry(x), preserveTopology, dTolerance):
st_simplify does not correctly simplify longitude/latitude data, dTolerance
needs to be in decimal degrees
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
sf_use_s2(TRUE)Spherical geometry (s2) switched on
Read Point Count
point_count <- read_csv(here::here('data', 'point_count.csv')) %>%
clean_names()New names:
Rows: 683669 Columns: 39
── Column specification
──────────────────────────────────────────────────────── Delimiter: "," chr
(19): GlobalUniqueIdentifier, ProjectCode, ProjectName, LocalityID, Stu... dbl
(16): ...1, SamplingUnitId, ParentSamplingUnitId, DecimalLatitude, Deci... lgl
(2): ProportionAreaSurveyed, InFocalHabitat date (1): ObservationDate time (1):
Time
ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
Specify the column types or set `show_col_types = FALSE` to quiet this message.
• `` -> `...1`
Read area search
area_search <- read_csv(here::here('data', 'area_search.csv')) %>%
clean_names()Rows: 163889 Columns: 37
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (15): GlobalUniqueIdentifier, ProjectCode, ProjectName, LocalityID, Stu...
dbl (15): SamplingUnitId, ParentSamplingUnitId, DecimalLatitude, DecimalLon...
lgl (5): AreaSurveyed, DistanceFromObserver, Salinity, Temperature, DepthMin
date (1): ObservationDate
time (1): Time
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Join both and make Geodatabase
point_area_geo <- full_join(area_search, point_count) %>%
st_as_sf(coords = c("decimal_longitude", "decimal_latitude"), crs = 4326) Joining with `by = join_by(global_unique_identifier, project_code,
project_name, locality_id, study_area, sampling_unit_id,
parent_sampling_unit_id, sampling_unit_type_name, decimal_latitude,
decimal_longitude, visit, protocol_code, observation_date, year_collected,
month_collected, day_collected, julian_day, julian_day_ve, time, collector,
scientific_name, common_name, species_code, phylogen_order, taxon_id,
distance_from_observer, observation_count, no_observations,
record_permissions)`
Leaflet layers
# Create color palettes
gap_pal <- colorFactor(
palette = c("green", "yellow", "orange", "red"), # Adjust colors as needed
domain = gap_clean$gap_sts # Or whatever your column name is
)
# Map with colored polygons
leaflet_map <- leaflet() %>%
addProviderTiles(providers$Esri.WorldTerrain, group = "ESRI Terrain") %>%
addMiniMap(toggleDisplay = TRUE) %>%
# GAP Status with colors
addPolygons(data = gap_clean,
group = 'GAP Status',
fillColor = ~gap_pal(gap_sts), # Replace gap_sts with your column
color = "black",
weight = 1,
fillOpacity = 0.6) %>%
# Vegetation Type with colors
addRasterImage(veg_type,
group = 'Vegetation Type',
opacity = 0.8) %>%
# Add legends
addLegend(pal = gap_pal,
values = gap_clean$gap_sts, # Replace with your column
title = "GAP Status",
position = "bottomleft",
group = 'GAP Status') %>%
addLayersControl(
overlayGroups = c('GAP Status', 'Vegetation Type'),
options = layersControlOptions(collapsed = TRUE)
) %>%
leaflet.extras::addResetMapButton()
leaflet_map